#ifndef EMBEDEDGRID_H
#define EMBEDEDGRID_H

#include <vector>
#include <sstream>
#include "Core/VecCompare.h"
#include "Core/Tensor.h"
#include "FiniteDiff/RectDomain.h"

template <int Dim>
class EmbeddedGrid : public RectDomain<Dim>
{
public:
  using BaseClass = RectDomain<Dim>;
  using rVec      = Vec<Real,Dim>;
  using iVec      = Vec<int,Dim>;
  template <class T> using vector = std::vector<T>;

  enum NodeType { Exterior = -1, Bdry = 0, Interior = 1 };

  template <class T_GeoInfo>
  EmbeddedGrid(const BaseClass& aRd,
               const T_GeoInfo &geoInfo,
               bool aUseOuterRect,
               Real tinyDist,
               Real aTol);

  // =================================================
  EmbeddedGrid(const EmbeddedGrid &)            = delete;
  EmbeddedGrid& operator=(const EmbeddedGrid &) = delete;
  EmbeddedGrid(EmbeddedGrid &&)                 = delete;
  EmbeddedGrid& operator=(EmbeddedGrid &&)      = delete;

  void dump(std::ostream& os) const;

public:
  using BaseClass::lo;
  using BaseClass::hi;
  using BaseClass::getDelta;
  using BaseClass::spacing;
  using BaseClass::getGhostedBox;
  using BaseClass::getNumGhost;

public:
  // geometry info
  Tensor<int,Dim>    nodeLabels;
  Tensor<rVec,Dim>   bdryPoints;
  Tensor<rVec,Dim>   outNormals;
  // boundary nodes and ghost nodes
  Tensor<iVec,Dim>   associatedGhost;
  vector<iVec>       bdryNodes;
  vector<iVec>       ghostNodes;
  // for direct solvers
  Tensor<int,Dim>    unknownIndices;
  int                numUnknowns;
  // others
  bool               useOuterRect;
protected:
  Box<Dim> valid;
  Box<Dim> gbx;
};

//============================================================

template <int Dim>
template <class T_GeoInfo>
inline
EmbeddedGrid<Dim>::EmbeddedGrid(const BaseClass &aRd,
                                const T_GeoInfo &geoInfo,
                                bool aUseOuterRect,
                                Real tinyDist,
                                Real aTol)
    : BaseClass(aRd),
      useOuterRect(aUseOuterRect)
{
  valid = *this;
  gbx = getGhostedBox();
  nodeLabels.resize(valid);
  bdryPoints.resize(valid);
  outNormals.resize(valid);
  associatedGhost.resize(valid);
  unknownIndices.resize(valid);

  // 1. Mark the interior/exterior nodes.
  tinyDist *= norm(spacing(), 2);
  ddfor(valid, [&](iVec idx) {
    rVec q = (idx + getDelta()) * spacing();
    rVec closest, normal;
    bool contained;
    geoInfo(q, closest, normal, contained, aTol);
    Real dist = norm(closest - q, 2);
    if(dist < tinyDist)
      nodeLabels(idx) = NodeType::Exterior;
    else
      nodeLabels(idx) = (contained) ? (NodeType::Interior) : (NodeType::Exterior);
    bdryPoints(idx) = closest;
    outNormals(idx) = normal;
  });

  // 2. Mark the boundary nodes
  numUnknowns = 0;
  ddfor(valid, [&](iVec idx) {
    if(nodeLabels(idx) == NodeType::Exterior)
      return;
    unknownIndices(idx) = numUnknowns++;
    // Check if this should be a boundary node
    for(int d = 0; d < Dim; ++d) {
      for(int side : {-1, 1}) {
        auto jdx = idx + iVec::unit(d) * side;
        if (valid.contain(jdx) && nodeLabels(jdx) == NodeType::Exterior) {
          nodeLabels(idx) = NodeType::Bdry;
          break;
        }
      }
      if(nodeLabels(idx) == NodeType::Bdry) {
        bdryNodes.push_back(idx);
        break;
      }
    }
  }, false);

  // 3. Associate each boundary node with a ghost node
  Tensor<int, Dim> ghostFlag(valid);
  ghostFlag = 0;
  for(const auto &idx : bdryNodes) {
    bool foundGhost = false;
    auto n = outNormals(idx);
    for(int k = 1; k <= 4; ++k) {
      iVec jdx = static_cast<iVec>(idx + getDelta() + n * k);
      if(valid.contain(jdx) && nodeLabels(jdx) == NodeType::Exterior && ghostFlag(jdx) == 0) {
        foundGhost = true;
        ghostFlag(jdx) = 1;
        associatedGhost(idx) = jdx;
        ghostNodes.push_back(jdx);
        break;
      }
    }
    if(!foundGhost) {
      std::ostringstream oss;
      oss << "Unable to assign a ghost node to " << idx << ". ";
      throw std::runtime_error(oss.str());
    }
  }
}

#endif // EMBEDEDGRID_H
